/*
MEX version of the gridding loop for the 2D NUFFT

Optimized version of nufftdauxmx_ref.cpp.
The iplemented optimizations are:
1) Move the varaible k to the outermost loop.
2) Remove unnecessary access to memory. For example, access mu only from the outer loop.
3) The arrays Px and Py are transposed compared to nufftdauxmx_ref, to access memory in columns.
4) The loops for j1,j2 now go from 0 to q instead of -q/2 to q/2.
5) Instad of computing the indices to Px, Py each time, we maintain the indices idxPx,idxPy, and just update them each iteration.

Example:
	tau=nufft3dauxmx(n,M,m,q,mu,Px.',Py.',alpha);
Note the transpose of Px and Py

Yoel Shkolnisky, January 2010.
*/

#include <mex.h>
#include <math.h>
#include <string.h>

//multithreading is based on code taken from here:
//http://www.mathworks.com/matlabcentral/fileexchange/21233-mex-multithreading
//http://www.yolinux.com/TUTORIALS/LinuxTutorialPosixThreads.html

#if defined __linux__ || defined __APPLE__
	#include <pthread.h>
#elif _WIN32 //defined for both 32 and 64 bit windows platform
	// undef needed for LCC compiler
	#undef EXTERN_C
	#include <windows.h>
	#include <process.h> 
#endif

#include "MEXUtilities.h"

#define for if(0);else for

typedef struct Data
{
	int threadN;
	int n;
	int M;
	int m;
	int q;
	double * mu;
	double * Px;
	double * Py;
	double * alpha_pr;
	double * alpha_pi;
	double * tau_pr;
	double * tau_pi;
	int ThreadID;
} DataType;

// The function which is multi threaded. 
// When ThreadID is k (1...feature('numCores')) it calculates the output of (k-1)/n ... k/n
// fraction out of the n required iterations
#if defined __linux__ || defined __APPLE__
	void * calculate(void * Args_in)
#elif _WIN32 //defined for both 32 and 64 bit windows platform
    unsigned __stdcall calculate(void * Args_in)
#endif
{
	DataType *Args = static_cast<DataType*>(Args_in);
    // Get all needed function variables from the double pointer array
 	int threadN,n,M,m,q,ThreadID;
	double *mu,*Px,*Py,*alpha_pr,*alpha_pi;	
	double *tau_pr, *tau_pi; // output variable
	mxComplexity complexFlag; // Is the input array alpha complex?
	threadN = Args->threadN;
	n = Args->n;
    M = Args->M;
    m = Args->m;
    q = Args->q;
	mu = Args->mu;
    Px = Args->Px;
    Py = Args->Py;
	alpha_pr = Args->alpha_pr;
    alpha_pi = Args->alpha_pi;
	complexFlag=mxREAL;
	if (alpha_pi) {
		complexFlag=mxCOMPLEX;
	}
    tau_pr = Args->tau_pr;
    tau_pi = Args->tau_pi;
    ThreadID = Args->ThreadID;

	long k,j1,j2;
	// We regrid each point omega(k) to (q+1)^2 regular points around it.
	// The regular grid points are the q+1 points closest to round(m*omega(k)).
	// The variables startinggridpointx, startinggridpointy hold the coordinates of first 
	// grid point of the regular grid around omega(k), in the x and y directions, respectively.
	// The coordinates are given in terms of the padded volume tau of size (m*M)^2.
	long startinggridpointx,startinggridpointy; 
	// Current grid point processed. "gp" stands for "grid point". (gpx,gpy) is the current 
	// grid point of the (q+1)^2 regular square around omega(k), at which we compute the conribution 
	// of the point omega(k).
	long gpx,gpy;
	// Weights to be used at the point (gpx,gpy). idxPx is the weight used to compute the contibution
	// of the point omega(k) at the point gpx, idxPy has a similar role.
	long idxPx,idxPy;	
	long gpxm,gpym; // The points gpx,gpy modulu m*M.
	long idx; // The linear index of the point (gpxm,gpym) in an array of size (m*M)^2.
	// Absolute value of the lowest index in a zero-centered array. For example, "2" if we have a 5 elements 
	// array, since the indices are -2,1,0,1,2.  This variable is used to translate indices of zero-centered 
	// arrays to indices of zero-based arrays.
	long offset; 
	double w,w1; // Temporary variables for computing the gridding weights.
	long mM; // To avoid repeated computations of m*M.
	double akr,aki; // Aviod accessing alpha_pr and alpha_pi for each j1,j2. Access it only once.

	// Gridding loop.
	mM=(long)(m*M);
	offset=(long)ceil((mM-1.0)/2.0);

	for (k=ThreadID; k<n; k=k+threadN) {
		startinggridpointx=(long)mu[k]-q/2;
		startinggridpointy=(long)mu[n+k]-q/2;
		gpx=startinggridpointx+offset;
		idxPx=k*(q+1);
		akr=alpha_pr[k];
		if (complexFlag) {
			aki=alpha_pi[k];
		}

		for (j1=0; j1<=q; j1++) {
			gpxm=matlab_mod(gpx,mM);
			w1=Px[idxPx];
			gpx++;
			idxPx++;
			gpy=startinggridpointy+offset;
			idxPy=k*(q+1);

			for (j2=0; j2<=q; j2++) {
				gpym=matlab_mod(gpy,mM);
				idx=gpym*mM+gpxm;
				w=w1*Py[idxPy];
				tau_pr[idx]+=w*akr;
				if (complexFlag) {
					tau_pi[idx]+=w*aki;
				}
				gpy++;
				idxPy++;
			}
		}
	}
	
	#ifdef _WIN32 //defined for both 32 and 64 bit windows platform
        // explicit end thread, helps to ensure proper recovery of resources allocated for the thread
        _endthreadex( 0 );
    #endif
    return 0;
}

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
	// Parameters: threadN,n,M,m,q,mu,Px,Py,alpha
	// Returns: tau
	int threadN,n,M,m,q;
	long mM;
	double *mu,*Px,*Py,*alpha_pr,*alpha_pi;	
	double **tau_pr, **tau_pi; // output variable
	mxComplexity complexFlag; // Is the input array alpha complex?
	int dims[2];
	if (nrhs!=9) {
		mexErrMsgTxt("10 parameters are required: threadN,n,M,m,q,mu,Px,Py,alpha");
	}
	// Get Input parameters
    GetInteger(prhs[0],&threadN);
	GetInteger(prhs[1],&n);
	GetInteger(prhs[2],&M);
	GetInteger(prhs[3],&m);
	GetInteger(prhs[4],&q);
	mu=mxGetPr(prhs[5]);
	Px=mxGetPr(prhs[6]);
	Py=mxGetPr(prhs[7]);
	alpha_pr=mxGetPr(prhs[8]);
	complexFlag=mxREAL;
	alpha_pi=NULL;
	if (mxIsComplex(prhs[8])) {
		alpha_pi=mxGetPi(prhs[8]);
		complexFlag=mxCOMPLEX;
	}
	mM=(long)(m*M);
	// Allocate output variable
	dims[0]=m*M; dims[1]=m*M;
	mxArray ** tau = (mxArray **)mxMalloc(threadN * sizeof( mxArray * ));
	tau_pr = (double **)mxMalloc((threadN+1)* sizeof( double * ));
	tau_pi = (double **)mxMalloc((threadN+1)* sizeof( double * ));
    for (int i = 0 ; i < threadN + 1 ; i++) {
		if (i == threadN)
		{
			plhs[0]=mxCreateNumericArray(2,dims,mxDOUBLE_CLASS,complexFlag);
			tau_pr[i]=mxGetPr(plhs[0]);
			tau_pi[i]=mxGetPi(plhs[0]);
		} else {
			tau[i]=mxCreateNumericArray(2,dims,mxDOUBLE_CLASS,complexFlag);
			tau_pr[i]=mxGetPr(tau[i]);
			tau_pi[i]=mxGetPi(tau[i]);
		}
    }

	// store all needed function variables
    DataType * ThreadArgs = (DataType *) mxMalloc(threadN * sizeof( DataType )); 
    #if defined __linux__ || defined __APPLE__
		pthread_t * ThreadList = (pthread_t *) mxMalloc(threadN * sizeof( pthread_t ));
    #elif _WIN32 //defined for both 32 and 64 bit windows platform
		// Handles to the worker threads
		// Reserve room for handles of threads in ThreadList
		HANDLE * ThreadList = (HANDLE *) mxMalloc(threadN * sizeof( HANDLE ));
	#endif
	for (int i = 0 ; i < threadN ; i++) {
		// create separate thread and start the calculate() function.
		ThreadArgs[i].threadN=threadN;
		ThreadArgs[i].n=n;
		ThreadArgs[i].M=M;
		ThreadArgs[i].m=m;
		ThreadArgs[i].q=q;
		ThreadArgs[i].mu=mu;
		ThreadArgs[i].Px=Px;
		ThreadArgs[i].Py=Py;
		ThreadArgs[i].alpha_pr=alpha_pr;
		ThreadArgs[i].alpha_pi=alpha_pi;
		ThreadArgs[i].tau_pr=tau_pr[i];
		ThreadArgs[i].tau_pi=tau_pi[i];
		ThreadArgs[i].ThreadID=i;
		#if defined __linux__ || defined __APPLE__
			pthread_create( &(ThreadList[i]), NULL, &calculate, (void *)&(ThreadArgs[i]));
		#elif _WIN32 //defined for both 32 and 64 bit windows platform
			ThreadList[i] = (HANDLE)_beginthreadex( NULL, 0, &calculate, (void *)&(ThreadArgs[i]) , 0, NULL );
		#endif
	}
    #if defined __linux__ || defined __APPLE__
		for (int i = 0 ; i < threadN ; i++) {
			pthread_join( ThreadList[i], NULL);
		}
    #elif _WIN32 //defined for both 32 and 64 bit windows platform
        // Normally, the next code line is used to wait for threads to finish
        WaitForMultipleObjects(threadN, ThreadList, true, INFINITE);
    #endif
	for (int j = 0 ; j < threadN ; j++) {
		for (int i = 0 ; i < dims[0]*dims[1] ; i++)
		{
			tau_pr[threadN][i]=tau_pr[threadN][i]+tau_pr[j][i];
		}
	}
	if (complexFlag==mxCOMPLEX) {
		for (int j = 0 ; j < threadN ; j++) {
			for (int i = 0 ; i < dims[0]*dims[1] ; i++)
			{
				tau_pi[threadN][i]=tau_pi[threadN][i]+tau_pi[j][i];
			}
		}
	}
    #ifdef _WIN32 //defined for both 32 and 64 bit windows platform
        // Explicit destroy the Mutex and Thread objects, helps to ensure proper recovery of resources.
		for (int i = 0 ; i < threadN ; i++) {
			CloseHandle( ThreadList[i] );
		}
    #endif
}

// Matlab reference code for the above main loop:
//for j1=-q/2:q/2
//    idx1=mu(:,1)+j1;
//    idx1=mod(idx1+offset1,m*M);
//    w1=Px(:,j1+q/2+1);
//    for j2=-q/2:q/2
//        idx2=mu(:,2)+j2;
//        idx2=mod(idx2+offset1,m*M);
//        w2=Py(:,j2+q/2+1);
//        w12=w1.*w2;
//        for j3=-q/2:q/2
//            idx3=mu(:,3)+j3;
//            idx3=mod(idx3+offset1,m*M);
//            idx=idx3.*(m*M)^2+idx2.*m*M+idx1+1; % convert to linear index
//            tau(idx(:))=tau(idx(:))+...
//                w12.*Pz(:,j3+q/2+1).*alpha(:);
//        end
//    end
//end
